package com.bodystm.algorithm;
/**
 * 血压血氧低通滤波
 * @author ehl
 *	cMax cCOMPLEX为内部定义的两个类（对应C++中的结构体）
 */
/*****************************************************************
File name:	functions for fast fourier transform under low perfusion condition
Author:	Siyang Liang
Version:1.0.0.20170329
Date: 2017/03/29
Description: low perfusion fft algorithm functions are provided
1. Date:2017/03/29
Author: ID: Siyang Liang
2. ...
*****************************************************************/
public class LowPerf {
	/**
	 * 两个算法都是需要1024个点计算一次，之后将后512个点保存下来，除第一次外后续的每一次都要用当上一次的最后512个点。
	 */
	public static final int SPO2CALCBUFFLEN=512;//1024;
	public int curNum=0;//当前累积的值的数量
	public float[] m_dIrf_AC=new float[SPO2CALCBUFFLEN];
	public float[] m_dIrf_DC=new float[SPO2CALCBUFFLEN];
	public float[] m_dRed_AC=new float[SPO2CALCBUFFLEN];
	public float[] m_dRed_DC=new float[SPO2CALCBUFFLEN];
	public int pushValue(double[] arr,double value){
		return -1;
	}
	/*public void resetDataArr(){
		for(int i = 0; i < 512; i++)
        {
            m_dRed_AC[i] = m_dRed_AC[i+512];
            m_dRed_DC[i] = m_dRed_DC[i+512];
            m_dIrf_AC[i] = m_dIrf_AC[i+512];
            m_dIrf_DC[i] = m_dIrf_DC[i+512];
        }
		curNum = 512;
	}*/
	public void resetDataArr(){
		curNum = 0;
	}
	/**
	* @describe: real time butterworth filtering, fc = 5Hz, order 5;
	* @param: input data point xn
	* @return: filtered datapoing yn
	* @author: Siyang Liang
	* @date: 2017/03/29
	*/
	double b[] = { 2.139615203814e-05, 1.069807601907e-04, 2.139615203814e-04, 2.139615203814e-04, 1.069807601907e-04, 2.139615203814e-05 };
	double a[] = { -4.187300047864, 7.069722752792, -6.009958148187, 2.570429302524, -4.422091823996e-01 };
	int CHANID=2;//红外和热红外||血压两种
	double bx[][] = new double[CHANID][6];
	double ax[][] =new double [CHANID][5];
	public double butter(double xn,int ChanID)
	{
		double atemp = 0,btemp = 0,yn = 0;
		bx[ChanID][5] = bx[ChanID][4];
		bx[ChanID][4] = bx[ChanID][3];
		bx[ChanID][3] = bx[ChanID][2];
		bx[ChanID][2] = bx[ChanID][1];
		bx[ChanID][1] = bx[ChanID][0]; //right shift b by 1;
		bx[ChanID][0] = xn;

		yn = (bx[ChanID][0] * b[0] + bx[ChanID][1] * b[1] + bx[ChanID][2] * b[2] + bx[ChanID][3] * b[3] + bx[ChanID][4] * b[4] + bx[ChanID][5] * b[5]) - (ax[ChanID][0] * a[0] + ax[ChanID][1] * a[1] + ax[ChanID][2] * a[2] + ax[ChanID][3] * a[3] + ax[ChanID][4] * a[4]);

		ax[ChanID][4] = ax[ChanID][3];
		ax[ChanID][3] = ax[ChanID][2];
		ax[ChanID][2] = ax[ChanID][1];
		ax[ChanID][1] = ax[ChanID][0]; //right shift a by 1
		ax[ChanID][0] = yn;
		return yn;
	}
	cCOMPLEX RedACComplex[]=new cCOMPLEX[1024];
	cCOMPLEX IredACComplex[]=new cCOMPLEX[1024];
	cCOMPLEX IredCorr[]=new cCOMPLEX[1024];
	cCOMPLEX BPMIred[]=new cCOMPLEX[2048];
	public LowPerf(){
		for (int i = 0; i < 1024; i++)	// initialize complex array for both light sources i.e. {real1,0,real2,0}
		{
	    	RedACComplex[i]=new cCOMPLEX();
	    	IredACComplex[i]=new cCOMPLEX();
	    	IredCorr[i]=new cCOMPLEX();
	    	BPMIred[i]=new cCOMPLEX();
	    	BPMIred[i+1024]=new cCOMPLEX();
		}
	}
	/**
	* @describe: new blood oxygen calculation algorithm for low perfusion condition
	* @param: 1024 length double vector 
	* @return: Object[]{SpO2, BPM,fPI}  int int float
	* @author: Siyang Liang
	* @date: 2017/03/29
	*/
	public Object[] LPSpO2Calc(double RedAC[], double RedDC[], double IredAC[], double IredDC[])
	{
		int SpO2; int BPM; float fPI;
		int i, N = 1024,Fbase = 4;
		double RDC = 0, IRDC = 0;
	    double tempSD = 0,diff,SD;
		float RAC, IRAC,R;
		float LocalMax,NextMax;
//		float P[3] = { -10.9329, -20.3176, 112.4622 };
	    float P[] = new float[]{ -18.7169f, -11.6142f, 110.4818f };
	    /*cCOMPLEX RedACComplex[]=new cCOMPLEX[1024];
		cCOMPLEX IredACComplex[]=new cCOMPLEX[1024];
		cCOMPLEX IredCorr[]=new cCOMPLEX[1024];
		cCOMPLEX BPMIred[]=new cCOMPLEX[2048];*/

	    for (i = 0; i < N; i++)	// initialize complex array for both light sources i.e. {real1,0,real2,0}
		{
	    	/*RedACComplex[i]=new cCOMPLEX();
	    	IredACComplex[i]=new cCOMPLEX();
	    	IredCorr[i]=new cCOMPLEX();
	    	BPMIred[i]=new cCOMPLEX();
	    	BPMIred[i+1024]=new cCOMPLEX();*/
	        RedACComplex[i].real = (float) RedAC[i];
	        RedACComplex[i].imag = 0;
	        IredACComplex[i].real = (float) IredAC[i];
	        IredACComplex[i].imag = 0;
	        RDC += RedDC[i];
	        IRDC += IredDC[i];
		}
		//////////////////////SpO2 Calculation Algorithm//////////////////////////
		RDC = RDC / N;
		IRDC = IRDC / N;
		fft(N, RedACComplex); //get the fft of both light signals
		fft(N, IredACComplex);
		ConjugateComplex(N, IredACComplex, IredCorr); //Auto correlation algorithm
		// finding the largest energy location of the auto correlation array
		LocalMax = IredCorr[4].real*IredCorr[4].real + IredCorr[4].imag*IredCorr[4].imag;
		for (i = 5; i <36 ; i++) //4 and 35 are the coresponding 0.5Hz and 4.25Hz location in the 1024 FFT resolution
		{	
			NextMax = IredCorr[i].real*IredCorr[i].real + IredCorr[i].imag*IredCorr[i].imag; // abs value of a complex number without sqrt
			if (LocalMax < NextMax)
			{
				LocalMax = NextMax;
				Fbase = i; // base frequency location
			}
		}
		RAC = RedACComplex[Fbase].real*RedACComplex[Fbase].real + RedACComplex[Fbase].imag*RedACComplex[Fbase].imag;
		RAC = (float) Math.sqrt(RAC);
		IRAC = IredACComplex[Fbase].real*IredACComplex[Fbase].real + IredACComplex[Fbase].imag*IredACComplex[Fbase].imag;
		IRAC = (float) Math.sqrt(IRAC);
		R = (float) ((RAC / RDC) / (IRAC / IRDC)); 
		SpO2 = (int)((R*R*P[0] + R*P[1] + P[2]) + 0.5); //return of theSpO2 value
//	    *fPI = (IRAC / IRDC);
	    for(i=0;i<N-1;i++)
	    {
	         diff = 0 -IredAC[i];
	         if(diff < 0 )
	         {
	             diff = -diff;
	         }
	         diff = diff*diff;
	         tempSD += diff;
	    }
	    SD = Math.sqrt(tempSD/N)*3.0;
	    fPI = (float) ((SD / IRDC)*100);

		/////////////////////////BPM Calculation Algorithm////////////////////////
	    //ifft(N, IredACComplex); //inverse fft
	    for (i = 0; i < 2048; i++)	//complext number initialization
		{
	        if(i < N)
	        {
	            BPMIred[i].real = (float) IredAC[i];
	        }
	        else
	        {
	            BPMIred[i].real = 0;
	        }
	        BPMIred[i].imag = 0;
	    }
	    fft(2048, BPMIred);
		LocalMax = BPMIred[8].real*BPMIred[8].real + BPMIred[8].imag*BPMIred[8].imag;

		for (i = 9; i <71; i++) //same algorithm as the autocorrelation, to find the largest energy of the BPM with higher resolution
		{
			NextMax = BPMIred[i].real*BPMIred[i].real + BPMIred[i].imag*BPMIred[i].imag; // abs value of a complex number without sqrt
			if (LocalMax < NextMax)
			{
				LocalMax = NextMax;
				Fbase = i; // base frequency location
			}
		}
	    BPM = (int)((Fbase-1)*0.060738 * 60 + 0.5); // return of the BPM value
	    return new Object[]{SpO2, BPM,fPI};
	}


	/**
	* @describe: Slop Sum Function (SSF) Calculation
	* @param: input data point segment, threshold
	* @return: none
	* @author: Siyang Liang
	* @date: 2017/04/19
	*/

	double SSFBuffer[] = new double[16];
	double Threshold = -1;
	int BPMPre = -1;
	public int SSF_Calc(double IredData[])
	{
	    double SSFTemp = 0,Temp;
	    double temp = 0, NewThreshold;
	    int i,j;
	    int BPM;
	    double SSF[]=new double[1024];

	    for (i = 0; i < 1024; i++)
	    {
	        //left shift
	        for (j = 0; j < 15; j++)
	        {
	            SSFBuffer[j] = SSFBuffer[j + 1];
	        }
	        SSFBuffer[15] = IredData[i];
	        for (j = 0; j < 15; j++)
	        {
	            Temp = -(SSFBuffer[j + 1] - SSFBuffer[j]);
	            if (Temp <= 0)
	            {
	                Temp = 0;
	            }
	            SSFTemp += Temp;
	        }
	        SSF[i] = SSFTemp;
	        SSFTemp = 0;
	    }

	    if (Threshold == -1)
	    {
	        for (i = 0; i < 1024; i++)
	        {
	            temp += SSF[i];
	        }
	        Threshold = (temp / 1024)*2.5;
	    }
	    Object[] result=SSF_BPM(SSF, Threshold);
	    BPM=(int) result[0];NewThreshold=(double) result[1];
	    if(BPM == 0 && BPMPre != -1)
	    {
	        BPM = BPMPre;
	        NewThreshold = -1;
	    }
	    BPMPre = BPM;
	    Threshold = NewThreshold;
	    return BPM;
	}

	/**
	* @describe: Detect pulse wave based on SSF detection
	* @param: input SSF result, threshold, pointer BPM, pointer NewThreshold
	* @return: none
	* @author: Siyang Liang
	* @date: 2017/04/19
	*/
	//[BPM,NewThreshold]
	private Object[] SSF_BPM(double SSF[], double Threshold)
	{
		int BPM;double NewThreshold;
	    //double SegMax = -1;
	    double SSFMax;
	    double BPMSum = 0;
	    int XLocTemp = -1;
	    int k = 0, i = 0;
	    int XLoc[] = new int[40];
	    int NumLoc,flag = 0;

	    while (k < 1023)
	    {
	        if (SSF[k] > Threshold &&  SSF[k + 1] - SSF[k] > 0)
	        {
	            SSFMax = SSF[k];
	            while (SSF[k] > Threshold)
	            {
	                if (SSFMax < SSF[k + 1])
	                {
	                    SSFMax = SSF[k + 1];
	                    XLocTemp = k + 1;
	                }
	                k++;
	                if (k>=1023)
	                {
	                    if (SSF[k] - SSF[k - 1] > 0) // if epoch ends before reaching the local max point, flag is 1, and will exclude at final calculation
	                    {
	                        flag = 1;
	                    }
	                    break;
	                }

	            }
	            /*if (SegMax < SSFMax)
	            {
	            SegMax = SSFMax;
	            }*/
	            XLoc[i] = XLocTemp;
	            Threshold = 0.6*SSFMax;
	            i++;
	        }
	        k++;

	    }
	    NumLoc = i - 1;
	    if (NumLoc < 1) // no peak detected
	    {
	        BPM = 0;
	        NewThreshold = -1;
	        return new Object[]{BPM,NewThreshold};
	    }


	    int count = 0;
	    for (i = 0; i < NumLoc-flag; i++) //if flag is 1, loop will ignore the last XLoc value since its not real local max
	    {
	        if (XLoc[i] != -1)
	        {
	            BPMSum += XLoc[i + 1] - XLoc[i];
	            count++;
	        }
	    }
	    if (BPMSum == 0)
	    {
	        BPM = 0;
	        NewThreshold = -1;
	        return new Object[]{BPM,NewThreshold};
	    }
	    BPM = (int) (60 / ((BPMSum / count)*0.008)+0.5);
	    NewThreshold = 0.6*Threshold;
	    return new Object[]{BPM,NewThreshold};
	}


	/**
	* @describe:
	* @param:
	* @return: 0:disable    1:enable
	* @author: Siyang Liang
	* @date: 2017/04/21
	*/
	int SPO2Pre = -1; 
	int[] MAFBuffer = new int[5];
	int Start = 1;
	int State_Detect(int SPO2)
	{
	    int diff,BufferLen = 5,tempSum = 0,State = 1;
	    int i,count;

	    if (SPO2Pre != -1)
	    {
	        diff = Math.abs(SPO2Pre - SPO2); //calculate difference between two calculated spo2 values
	        if (Start <= BufferLen) //set Buffer values to all diff value when the program initializes
	        {
	            for (i = 0; i < BufferLen; i++)
	            {
	                MAFBuffer[i] = diff;
	            }
	            Start++;

	        }

	        for (i = 0; i < BufferLen-1; i++)	//left shift
	        {
	            MAFBuffer[i] = MAFBuffer[i + 1];
	            tempSum += MAFBuffer[i+1];
	        }
	        MAFBuffer[BufferLen - 1] = diff;
	        tempSum += diff;
	        tempSum = tempSum / BufferLen;
	    }
	    if (tempSum > 3.5) //MAF threshold
	    {
	        State = 0;
	    }
	    else
	    {
	        State = 1;
	    }
	    SPO2Pre = SPO2;
	    return State;
	}
	public class cCOMPLEX{
		public float real;
		public float imag;
		/*public float getReal() {
			return real;
		}
		public void setReal(float real) {
			this.real = real;
		}
		public float getImag() {
			return imag;
		}
		public void setImag(float imag) {
			this.imag = imag;
		}*/
	}
	public class cMax{
		public int location;	//location of the max value in an array
		public float value;	//value of the max
	}
	//=============================BODataProcessFunc.cpp==============================
	/*****************************************************************
	File name:	functions for fast fourier transform
	Author:	Siyang Liang
	Version:1.0.0.20170204
	Date: 2017/02/04
	Description: functions needed for calculating fast fourier transform, complex number 
				is a struct defined in the header file which contains two elements, real 
				part (cCOMPLEX.real)and imaginary part (cCOMPLEX.imag).
	Others: FFT algorithm complete simulation works fine, real time test still needed.
			newly added BPM calculation algorithm, simulation complete, real time test needed.
	History: none
	1. Date:2017/02/07
	Author: ID: Siyang Liang
	Modification: changed the input of SpO2 calculation and BPM calculation from float to double
				  minor changes to ifmprove performance speed
				  adding the new BPM algorithm
				  **NOTE** further improvement in speed still possible
	2. ...
	*****************************************************************/
	private static final int CHAN_COUNT=2;
	double dBuffer[][]=new double[CHAN_COUNT][15];

	//////////////////////////FFT functions///////////////////////////
	/**
	* @describe: conjugate of a complex number	array
	* @param: length of the array, input complex number array, output complex number array
	* @return: none
	* @author: unknown
	* @date: 2017/02/04
	*/
	void ConjugateComplex(int n, cCOMPLEX in[], cCOMPLEX out[])
	{
		int i = 0;
		for (i = 0; i<n; i++)
		{
			out[i].imag = -in[i].imag;
			out[i].real = in[i].real;
		}
	}

	/**
	* @describe: absolute value of a complex number array
	* @param: complex number array, absolute value ouput array, input array length
	* @return: none
	* @author: unknown
	* @date: 2017/02/04
	*/
	void ComplexABS(cCOMPLEX f[], float out[], int n)
	{
		int i = 0;
		float t;
		for (i = 0; i<n; i++)
		{
			t = f[i].real * f[i].real + f[i].imag * f[i].imag;
			out[i] = (float) Math.sqrt(t);
		}
	}

	/**
	* @describe: addition of two complex numbers
	* @param: complex number a, complex number b, result complex number struct pointer c
	* @return: none
	* @author: unknown
	* @date: 2017/02/04
	*/
	void ComplexPlus(cCOMPLEX a, cCOMPLEX b, cCOMPLEX c)
	{
		c.real = a.real + b.real;
		c.imag = a.imag + b.imag;
	}

	/**
	* @describe: subtraction of two complex numbers
	* @param: complex number a, complex number b, result complex number struct pointer c
	* @return: none
	* @author: unknown
	* @date: 2017/02/04
	*/
	void ComplexSub(cCOMPLEX a, cCOMPLEX b, cCOMPLEX c)
	{
		c.real = a.real - b.real;
		c.imag = a.imag - b.imag;
	}

	/**
	* @describe: multiplication of two complex numbers 
	* @param: complex number a, complex number b, result complex number struct pointer c
	* @return: none
	* @author: unknown
	* @date: 2017/02/04
	*/
	void ComplexMulti(cCOMPLEX a, cCOMPLEX b, cCOMPLEX c)
	{
		c.real = a.real * b.real - a.imag * b.imag;
		c.imag = a.real * b.imag + a.imag * b.real;
	}

	/**
	* @describe: division of two complex numbers
	* @param: complex number a, complex number b, result complex number struct pointer c
	* @return: none
	* @author: unknown
	* @date: 2017/02/04
	*/
	void ComplexDiv(cCOMPLEX a, cCOMPLEX b, cCOMPLEX c)
	{
		c.real = (a.real * b.real + a.imag * b.imag) / (b.real * b.real + b.imag * b.imag);
		c.imag = (a.imag * b.real - a.real * b.imag) / (b.real * b.real + b.imag * b.imag);
	}

	private void SWAP(Object a,Object b){
		Object tempr=a;a=b;b=tempr; //ignore for now
	}

	/**
	* @describe: wn to Wn tranformation for recursive calculation
	* @param: number of iteration, iteration value, resulting complex number, flag indicator
	* @return: none
	* @author: unknown
	* @date: 2017/02/04
	*/
	void Wn_i(int n, int i, cCOMPLEX Wn, int flag)
	{
		Wn.real = (float) Math.cos(2 * Math.PI*i / n);
		if (flag == 1)
		{
			Wn.imag = (float) -Math.sin(2 * Math.PI*i / n);
		}
		else if (flag == 0)
		{
			Wn.imag = (float) -Math.sin(2 * Math.PI*i / n);
		}
	}

	/**
	* @describe: fast fourier transform (fft) code obtained from the internet
	* @param: N is the fft length, f[] is the data array
	* @return: none
	* @author:	unkonwn
	* @date: 2017/02/04
	*/

	/***************Note*****************/
	/*the input array cCOMPLEX f[] need   */
	/*to be a complext number struct with only the real part*/
	/*imaginary part should be 0 in this case */
	/*i.e. input signal = [ 1 2 3 4 5 8 ]*/
	/*input f[] should be cCOMPLEX f[6] = {1,0,2,0,3,0,4,0,5,0,8,0}*/
	void fft(int N, cCOMPLEX f[])
	{
		cCOMPLEX t=new cCOMPLEX(), wn=new cCOMPLEX();//ÖÐ¼ä±äÁ¿  
		int i, j, k, m, n, l, r, M;
		int la, lb, lc;
		/*----¼ÆËã·Ö½âµÄ¼¶ÊýM=log2(N)----*/
		for (i = N, M = 1; (i = i / 2) != 1; M++);
		/*----°´ÕÕµ¹Î»ÐòÖØÐÂÅÅÁÐÔ­ÐÅºÅ----*/
		for (i = 1, j = N / 2; i <= N - 2; i++)
		{
			if (i<j)
			{
				t = f[j];
				f[j] = f[i];
				f[i] = t;
			}
			k = N / 2;
			while (k <= j)
			{
				j = j - k;
				k = k / 2;
			}
			j = j + k;
		}

		/*----FFTËã·¨----*/
		for (m = 1; m <= M; m++)
		{
			la = (int) Math.pow((double) 2, m); //la=2^m´ú±íµÚm¼¶Ã¿¸ö·Ö×éËùº¬½ÚµãÊý       
			lb = la / 2;    //lb´ú±íµÚm¼¶Ã¿¸ö·Ö×éËùº¬µúÐÎµ¥ÔªÊý  
			//Í¬Ê±ËüÒ²±íÊ¾Ã¿¸öµúÐÎµ¥ÔªÉÏÏÂ½ÚµãÖ®¼äµÄ¾àÀë  
			/*----µúÐÎÔËËã----*/
			for (l = 1; l <= lb; l++)
			{
				r = (int) ((l - 1)*Math.pow((double) 2, M - m));
				for (n = l - 1; n<N - 1; n = n + la) //±éÀúÃ¿¸ö·Ö×é£¬·Ö×é×ÜÊýÎªN/la  
				{
					lc = n + lb;  //n,lc·Ö±ð´ú±íÒ»¸öµúÐÎµ¥ÔªµÄÉÏ¡¢ÏÂ½Úµã±àºÅ       
					Wn_i(N, r, wn, 1);//wn=Wnr  
					ComplexMulti(f[lc], wn, t);//t = f[lc] * wn¸´ÊýÔËËã  
					ComplexSub(f[n], t, f[lc]);//f[lc] = f[n] - f[lc] * Wnr  
					ComplexPlus(f[n], t, f[n]);//f[n] = f[n] + f[lc] * Wnr  
				}
			}
		}
	}

	/**
	* @describe: inverse fft
	* @param: ifft length	
	* @return: none
	* @author: unkonwn
	* @date: 2017/02/04
	*/
	void ifft(int N, cCOMPLEX f[])
	{
		int i = 0;
		ConjugateComplex(N, f, f);
		fft(N, f);
		ConjugateComplex(N, f, f);
		for (i = 0; i<N; i++)
		{
			f[i].imag = (f[i].imag) / N;
			f[i].real = (f[i].real) / N;
		}
	}
	//////////////////////////////Tool Functions////////////////////////////////
	/**
	* @describe: finding the max of an array
	* @param: array pointer, array length
	* @return: maxmum value
	* @author: Siyang Liang
	* @date: 2017/02/04
	*/
	float Max(float[] p,int startIndex,int N)
	{
		int i;
		float fLocalMax = p[startIndex];
		for (i = 1; i < N-1; i++)
		{
			fLocalMax=fLocalMax < p[startIndex+i] ? p[startIndex+i] : fLocalMax;
		}
		return fLocalMax;
	}

	/**
	* @describe: finding the max of an array
	* @param: array pointer, array length
	* @return: maximum struct containing of value and location of the max
			   i.e. Max.value and Max.location	
	* @author: Siyang Liang
	* @date: 2017/02/07
	*/
	cMax StructMax(float[] p,int startIndex, int N)
	{	
		cMax CMax=new cMax();
		int i;
		//float fLocalMax = *p;
		CMax.value = p[startIndex];
		for (i = 1; i < N - 1; i++)
		{
			if (CMax.value < p[startIndex+i])
			{
				CMax.value = p[startIndex+i];
				CMax.location = i;
			}
			else
			{
				continue;
			}
		}
		return CMax;
	}

	/**
	* @describe: sum of an array
	* @param: array pointer, array length
	* @return: result of sum
	* @author: Siyang Liang
	* @date: 2017/02/04
	*/
	int Sum(int[] pt,int startIndex, int len)
	{
		int result = 0;
		for (int i = 0; i < len; i++)
		{
			result += pt[startIndex+i];
		}
		return result;
	}

	/**
	* @describe: blood oxygen sturation
	* @param: float type red light signal array, float type infrared light array, NOTE: both input array lengths need to be 256
	* @return: blood oxygen saturation level in percent
	* @author: Siyang Liang
	* @date: 2017/02/04
	*/

	int BloodOxygenCalc_(double dRedSignal[], double dIredSignal[], float[] pi)
	{
	    cCOMPLEX CComplexRedSignal[]=new cCOMPLEX[256];	//struct element of complex array for fft calculation
	    cCOMPLEX CComplexIredSignal[]=new cCOMPLEX[256];

	    int i, N = 256;
	    float fFFTRedResult[]=new float[256], fFFTIredResult[]=new float[256];	//abusolute value of the resulting fft of the signal
	    float fRedAC, fIredAC;	//AC component of red and infrared light signal
	    float fRedDC, fIredDC;	//DC component of red and infrared light signal
	    float fR;	//R value
//	    double dP[3] = { -15.3072, -14.1771, 110.4211};	//stored regression model obtained previously
	    double dP[] = new double[]{ -18.4740, -11.9697, 110.5819};	//stored regression model obtained previously


	    for (i = 0; i < N; i++)	// initialize complex array for both light sources i.e. {real1,0,real2,0}
	    {
	        CComplexRedSignal[i].real = (float) dRedSignal[i];
	        CComplexRedSignal[i].imag = 0;
	        CComplexIredSignal[i].real = (float) dIredSignal[i];
	        CComplexIredSignal[i].imag = 0;
	    }
	    //FFT calculation and computing the abusolute value
	    fft(N, CComplexRedSignal);
	    ComplexABS(CComplexRedSignal, fFFTRedResult, N);
	    fft(N, CComplexIredSignal);
	    ComplexABS(CComplexIredSignal, fFFTIredResult, N);
	    //AC and DC components extraction and R value calculation based on algorithm
	    fRedAC = 2 * Max(fFFTRedResult,2, 16) / N;
	    fIredAC = 2 * Max(fFFTIredResult,2, 16) / N;
	    fRedDC = fFFTRedResult[0] / N;
	    fIredDC = fFFTIredResult[0] / N;
	    pi[0] = fIredAC / fIredDC * 400;
	    fR = (fRedAC / fRedDC) / (fIredAC / fIredDC);

	    return (int) ((fR*fR*dP[0] + fR*dP[1] + dP[2]) + 0.5);	//returning the final result
	}

	/**
	* @describe: beat per minute (BPM) calculation based on FFT
	* @param: input signal, it can be either red light signal or infrared light
			  Note: the input signal should be without baseline wander (after median filter or HPF)	
	* @return: BPM value    1/15, 1/15, 1/15, 1/15, 1/15, 1/15, 1/15, 1/15, 1/15, 1/15, 1/15, 1/15, 1/15, 1/15, 1/15
	* @author: Siyang Liang
	* @date: 2017/02/07
	*/
	int BPMCalc(double dSignal[])
	{
		cCOMPLEX CfSignal[]=new cCOMPLEX[2048];	//trailing the original data with padding zeros
		cMax CLocal;

	    //memset(CLocal, 0, sizeof(cMax));
	    //memset(CfSignal, 0 , sizeof(cCOMPLEX));

		int i;
		float fFFTABS[]=new float[2048];

	    for (i = 0; i < 2048; i++)	//padding signal with trailing zeros
	    {
	        if(i < 512)
	        {
	        CfSignal[i].real = (float) dSignal[i];
	        CfSignal[i].imag = 0;
	        }
	        else
	        {
	            CfSignal[i].real = 0;
	            CfSignal[i].imag = 0;
	        }
	    }

		//FFT algorithm for calculationg BPM
		fft(2048, CfSignal);//[0]
		ComplexABS(CfSignal, fFFTABS, 2048);
	    CLocal = StructMax(fFFTABS,7, 64);
	    return (int) ((CLocal.location + 7)*0.060135*60+0.5);
	}

	/**
	* @describe: weighted moving average filter (WMAF)
	* @param: input data point
	* @return: none
	* @author: Siyang Liang
	* @date: 2017/02/13
	*/
	double WMAF(double dInput, int nChanId)//, double *dOutput
	{
	    int i;
	    double dOutput = 0;

	    for (i = 0; i < 14; i++)
	    {
	        dBuffer[nChanId][i] = dBuffer[nChanId][i + 1];
	    }
	    dBuffer[nChanId][14] = dInput;
	    for (i = 0; i < 14; i++)
	    {
	        dOutput += dBuffer[nChanId][i];
	    }
	    dOutput = dOutput / 15;
	    return dOutput;
	}
}
